Ultra-strong spin–orbit coupling and topological moiré engineering in twisted ZrS2 bilayers

We predict that twisted bilayers of 1T-ZrS2 realize a novel and tunable platform to engineer two-dimensional topological quantum phases dominated by strong spin-orbit interactions. At small twist angles, ZrS2 heterostructures give rise to an emergent and twist-controlled moiré Kagome lattice, combining geometric frustration and strong spin-orbit coupling to give rise to a moiré quantum spin Hall insulator with highly controllable and nearly-dispersionless bands. We devise a generic pseudo-spin theory for group-IV transition metal dichalcogenides that relies on the two-component character of the valence band maximum of the 1T structure at Γ, and study the emergence of a robust quantum anomalous Hall phase as well as possible fractional Chern insulating states from strong Coulomb repulsion at fractional fillings of the topological moiré Kagome bands. Our results establish group-IV transition metal dichalcogenide bilayers as a novel moiré platform to realize strongly-correlated topological phases in a twist-tunable setting.

T wisted van der Waals heterostructures have recently emerged as an intriguing and highly tunable platform to realize unconventional electronic phases in two dimensions [1][2][3][4] . Spurred by the discovery of Mott insulation and superconductivity in twisted bilayer graphene 5,6 , remarkable progress in fabrication and twist-angle control has led to observations of correlated insulating states or superconductivity in a variety of materials, including trilayer and double-bilayer graphene, homo-and hetero-bilayers of twisted transition metal dichalcogenides (TMDs) [7][8][9][10][11][12][13][14][15][16][17][18] , and heterostructures at a twist on hexagonal boron nitride substrates 19,20 . At its heart, this rich phenomenology stems from electronic interference effects due to the moiré superlattice, which can selectively quench kinetic energy scales to realize almost dispersionless bands, permitting a twist angle controlled realization of regimes dominated by strong electronic interactions. At the same time, the drastic reduction of kinetic energy of the low-energy moiré bands implies straightforward gate-tunable access to a wide range of filling fractions, permitting wide-ranging experimental access to the phase diagrams of paradigmatic models of stronglycorrelated electrons 2 . Consequently, the putative realization of strongly-correlated electron physics in a tunable setting has garnered significant attention, resulting in growing experimental evidence for novel correlated phases, including unconventional superconductivity 4,21,22 .
Notably, and despite negligible intrinsic spin-orbit coupling in graphene, these were found to include topological states of matter. Here, the realization of the interaction-induced quantum anomalous Hall effect without external magnetic fields in twisted bilayer [23][24][25] and trilayer 26 graphene has spurred numerous proposals for more exotic fractionalized topological states of matter [27][28][29][30] , which however rely on a delicate interplay of spontaneous ferromagnetic order, valley polarization, and substrate engineering effects to induce the requisite nontrivial band topology. Generalizations to twisted transition-metal dichalcogenides have focused on telluride-based group VI compounds with 2H structure in the monolayer which exhibit an intrinsic quantum spin Hall effect 31 , with the quantum anomalous Hall effect recently observed 32 and similarly expected to emerge from spontaneous valley polarization 33,34 .
Central to the present work, we demonstrate for twisted bilayer ZrS 2 with 1T structure that the paradigm of twist-controlled suppression of the bare kinetic energy scales can be straightforwardly extended to instead promote spin-orbit coupling to constitute the dominant energy scale at low energies, opening up a new and exotic regime for experimental and theoretical investigation. Remarkably, we find that the two-component character of the valence band maximum in such two-dimensional group IV transition metal dichalcogenides enters in an essential manner, leading to the emergence of a clean moiré Kagome lattice with almost dispersionless quantum spin Hall bands at small twist angles. We demonstrate that this tunable realization of a ZrS 2 moiré heterostructure with strong spin-orbit coupling and strong interactions can therefore provide a robust and novel platform to probe the profound interplay of non-trivial band topology and electronic correlations, and shed light on elusive quantum phases beyond the purview of conventional condensed matter systems.

Results
Emergent Kagome moiré pattern in twisted ZrS 2 bilayers. ZrS 2 is a group IV transition metal dichalcogenide with an exfoliable layered structure. Different from group VI TMDs such as Mo 2 and WS 2 that normally adopt a 2H layered structure, ZrS 2 has a stable 1T structure 35 in its ground state as shown in Fig. 1a, without distorting into the 1T' structure 36,37 . Bulk ZrS 2 is a semiconductor with a band gap of 1.80 eV 38,39 and it remains semiconducting when thinned down to the monolayer 35,40 . In contrast to group-VI transition metal dichalcogenides such as MoS 2 with 2H structure, the valence band maximum in ZrS 2 and other group-IV transition metal dichalcogenides is located at Γ already in the monolayer and is composed of twofold degenerate chalcogen p x , p y orbitals. Spin-orbit coupling lifts their degeneracy and introduces a~100 meV gap [ Fig. 1b]. This property readily carries over to aligned bilayers with symmetric AA and AB stacking configurations [ Fig. 1c, d]; here, the valence band maximum at Γ follows from antibonding combinations of the out-of-plane chalcogen p x , p y orbitals. These are energetically separated from bonding combinations by~80−100 meV [ Fig. 1c, d], with a secondary local valence band maximum of p z orbitals furthermore located close to Γ and similarly detuned by~50 meV for AA stacking.
In twisted bilayers, the atomic interlayer registry interpolates continuously between local AA, AB, and BA alignment as a function of position and a moiré pattern with three-fold rotation symmetry forms [ Fig. 1e]. At sufficiently small twist angles, the energetic considerations for aligned bilayers discussed above immediately suggest that the top-most (highest energy) moiré valence bands should be similarly composed of antibonding p x , p y chalcogen orbitals. If spin-orbit coupling is neglected, these are degenerate at Γ in both AA and AB regions [ Fig. 1e] by virtue of rotation symmetry. However, the valence band edge differs between the two stackings, with the smooth interpolation between local alignments in the moiré unit cell encoded in an effective periodic scalar moiré potential V(r) [Fig. 1f]. Minima of V(r) is located at the AB and BA regions and form an effective honeycomb lattice. Notably, for the purposes of capturing the highest-energy moiré valence bands, the potential retains to an excellent approximation of the full sixfold rotation and mirror symmetries of the monolayer, even though the macroscopic crystal is chiral. This situation is in principle analogous to twisted bilayer MoS 2 41,42 , which hosts a series of almost dispersionless bands of Mo d z 2 orbital character on an emergent moiré honeycomb lattice.
Crucially, the loss of rotational symmetry away from local AA, AB stacking lifts the orbital degeneracy between p x , p y antibonding orbitals (in the absence of spin-orbit coupling), introducing a second energy scale into the problem. In stark contrast to 2H TMD bilayers, the two-component character of the Γ valley states enters in an essential manner. From symmetry considerations, their orbital splitting is expected to be maximal in three "domain wall" regions "X" per moiré unit cell in Fig. 1e, in which the transition-metal atoms of both layers form "stripes", with the rotational symmetry of the local stacking order reduced to C 2 . Contrary to the scalar moiré potential, maxima of orbital p x , p y splitting hence form a Kagome pattern [ Fig. 1g]. Remarkably, if the resulting energetic gain exceeds the scalar potential V(r), it becomes favorable for charge to migrate from the honeycomb AB/BA regions to "X" regions, realizing an emergent Kagome lattice of s-like moiré orbitals in a highly-tunable setting [ Fig. 1e].
Continuum model of twisted ZrS 2 . A minimal continuum model of this scenario readily follows from the above symmetry considerations asĤ whereĤ 0 describes the two-fold degenerate antibonding p x , p y chalcogen stateŝ ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-31604-w with the orbital degree of freedom represented via Pauli matriceŝ τ. Here, m ⋆ denotes the effective average band mass, and η ¼ m þ Àm À m þ þm À parametrizes the ratio of light (m − ) and heavy (m + ) hole p bands at Γ. Atomic spin-orbit interactionŝ lift the orbital degeneracy, opening up a gap at Γ as discussed in detail below. Here,σ z acts on spin. Central to the emergence of the Kagome lattice, the moiré potential acts nontrivially on the orbital pseudospin, and can generically be written as a Fourier expansion Here, n indexes the n-th moiré Brillouin zone. V n parameterizes the Fourier modes of the scalar potential in direct analogy to twisted WS 2 42 , with f ð0Þ n ðrÞ ¼ cosðb n;1 rÞ þ cosðb n;2 rÞ þ cosðb n;3 rÞ chosen to retain the full sixfold rotation symmetry and b n,i describing the three reciprocal lattice vectors i = 1, 2, 3 (related via C 3 rotations) to the nth Brillouin zone.
The pseudospinτ x ,τ z contributions to the potential are related in the presence of (approximate) mirror symmetry, with f ðxÞ n ðrÞ ¼ À ffiffi 3 p 2 cosðb n;1 rÞ þ ffiffi 3 p 2 cosðb n;3 rÞ and f ðzÞ n ðrÞ ¼ 1 2 cosðb n;1 rÞ À cosðb n;2 rÞ þ 1 2 cosðb n;3 rÞ. The salient physics is encoded already in the lowest harmonic-with b 1;1 ¼ ½2π; À2π= ffiffi ffi 3 ffiffi ffi 3 p =a 0 and a 0 the moiré lattice length, the scalar potential hosts two minima in the AB and BA regions at r ¼ ½1=2; ± 1=ð2 ffiffi ffi 3 p Þ. Conversely, the pseudospin potential that determines the splitting of p x , p y orbitals has three maxima in the Kagome X regions at r ¼ ½1=2; 0a 0 ; ½1=4; ffiffi ffi 3 p =4a 0 ; ½3=4; ffiffi ffi 3 p =4a 0 . Up to an overall energy scale, a minimal continuum model that includes only the first harmonic will therefore depend on just three dimensionless parameters η, 2Vm ? a 2 0 =_ 2 , 2V 0 m ? a 2 0 =_ 2 , where a 0 ∝ θ −1 scales with the twist angle. Local lattice relaxation effects are encoded in the higher harmonics of the potential; as the local stacking of AB and BA regions arise energetic favorable, lattice relaxation results in large domains with almost uniform AB or BA stacking [see domains highlighted with red and purple dashed lines in Fig. 1e] with the salient stacking variation in the X regions. These parameters can be obtained by fitting the band structure obtained from DFT calculations as described below. , and an energetically well-separated set of honeycomb bands with Dirac points at K, K 0 emerges at the top of the valence band. These retain a twofold orbital p x , p y character, with the degeneracy of the bands weakly broken due to orbital anisotropy η ≠ 0. This directly mirrors the low-energy band structure of twisted bilayer graphene, however with the two-orbital structure resulting from the p x , p y degeneracy of the constituent states at Γ as opposed to a valley degeneracy. The band structure of ZrS 2 monolayer hosts a valence band maximum at Γ composed of the p x , p y chalcogen orbital J = ± 3/2 states, with a spin-orbit splitting of~100 meV. c, d In aligned bilayers with AA or AB stacking, the valence band maximum at Γ (composed of of antibonding states of p x , p y symmetry) is energetically separated by~50meV from bonding p x , p y states and p z orbitals. Hence, antibonding p x , p y states compose the top-most moiré valence bands at small twist angles. e In small-angle twisted heterostructures of ZrS 2 , the interlayer alignment interpolates from AA to AB/BA stacking as a function of position in the moiré unit cell. The p x , p y antibonding states remain degenerate for both AA and AB/BA stacking in the absence of spin-orbit coupling. Their overall energetic shift is encoded in f as an effective scalar moiré superlattice potential, with minima that form a moiré honeycomb lattice. g Crucially, the loss of rotational symmetry away from local AA, AB stacking further splits the p x , p y states. This orbital splitting contribution to the moiré potential is maximal in three "domain wall" regions "X" in the moiré unit cell [(e)]; upon exceeding the scalar potential, electrons in the top-most moiré valence band form an emergent moiré Kagome lattice [dashed lines in e; guide to the eye].
However, already the next lower in energy (fifth) moiré valence band reveals upon closer inspection a charge density distribution with a Kagome pattern [ Fig. 2b, pattern (II)], localized in the "X" regions of the moiré unit cell [Fig. 1e]. These states gain energy from a finite pseudospin potential V 0 , which lifts them to higher energies: Beyond a critical V 0 , the fifth "Kagome" band and the bottom p x , p y honeycomb bands invert their energetic ordering at Γ. Consequently, the charge density distribution of the top p x , p y bands shifts from AB/BA honeycomb regions to "X" Ab initio characterization. The above behavior closely matches the results from large-scale ab initio calculations of the twisted moiré supercell, depicted in Fig. 2d, left column, for three representative twist angles [see Supplementary Note 1]. As the angle is reduced, a set of bands with a Kagome charge distribution at Γ splits off progressively from deeper valence bands. For the larger twist angles ≥2.28 ∘ that are still within computational reach for density functional calculations, this energetic separation is not yet sufficient to completely separate the Kagome bands of chalcogen antibonding p x , p y character from states with p z or bonding As the top-most valence states originate from p x , p y orbitals at small twist angles, spin-flip spin-orbit interactions are negligible and spin-z remains a good quantum number. Remarkably, this results in three almost dispersionless moiré bands that realize a novel Kagome topological quantum spin Hall insulator with spin Chern numbers C s ¼ ±1 for the first and third flat band [Fig. 2c, d]. In marked contrast to conventional topological materials, however, while superlattice interference quenches the kinetic energy scales, spin-orbit coupling λ soc enters as a bare atomic scale and hence becomes the dominant energy scale that governs the low-energy physics of the moiré valence bands in ZrS 2 . This highly-tunable materials realization of an "ultra-strong" spinorbit interaction regime in a moiré heterostructure constitutes a central result of this paper.
To model the emergent top-most flat topological moiré band in twisted ZrS 2 , we proceed with a fit of the pseudospin continuum theory [Eq. (1)] to the spin-orbit-coupled ab initio band structure for θ = 2.64 ∘ [Fig. 2d, middle-right panel]. As the minimal model of Eq. (1) does not account for bonding p x , p y , or p z states, the third-highest ab initio valence band (−50 meV below the valence band edge) is composed primarily of bonding p x , p y , and p z orbitals and is excluded from the fit. We note that this band separates energetically from the three Kagome moiré bands at lower twist angles. We obtain excellent agreement for the top two bands of p x , p y antibonding character using η = 0.33, m ⋆ = 0.27m 0 , λ soc = 57 meV, V 1 = 5.5 meV, V 0 1 ¼ À9:3 meV, V 2 = 11.5 meV, V 0 2 ¼ À5:1 meV. Scaling with twist angle similarly matches the ab initio band structure at 2.45 ∘ [Fig. 2d, bottom-right panel]. As expected, the top-most band is topologically nontrivial with spin Chern number C s ¼ ±1. Figure 2e compares the corresponding charge density distributions at Γ for ab initio and continuum model calculations; both exhibit comparable Kagome patterns as well as a competing band at lower energies with a ring-shaped charge pattern around the AA region, which similarly becomes energetically separated from Kagome bands at lower twist angles [ Fig. 2a].

Tight-binding description of emergent moiré Kagome bands.
A key advantage of the continuum theory is the possibility to study the behavior at small twist angles in a computationally feasible manner. Figure 3a depicts the bandwidth of the top-most topological moiré Kagome band, as well as the single-particle gap to the next deeper valence band, as a function of twist angle a 0~θ −1 . The bandwidth of the top-most topological band decreases exponentially with twist angle, whereas the ratio between bandwidth and band gap saturates below ≈ 2 ∘ and approaches one. Below this twist angle, the three Kagome bands become fully isolated in energy from deeper valence states [ Fig. 3g]. This immediately suggests a fruitful tight-binding parameterization at ultra-small angles, presuming that local lattice relaxation effects remain manageable. Results are shown in Fig. 3b for a tightbinding model depicted schematically in (c), but including up to 8th-neighbor hopping to ensure a good fit over all angles [see Supplementary Note 2]. For small angles ≪ 2 ∘ , the top three bands become well-captured by a nearest-neighbor Kagome tightbinding model with imaginary hoppings. Third-neighbor hopping t δ through the hexagons are leading corrections to this model and follow from the elliptical shapes of the charge density distribution at the Kagome "X" sites.
The sizable imaginary nearest-neighbor hopping [ Fig. 3b] is a direct consequence of the strong spin-orbit coupling limit and can be interpreted as a finite effective staggered magnetic flux through the elementary triangles of the Kagome lattice. It lifts the quadratic touching of flat and dispersive Kagome bands and opens up a gap at the Dirac points, realizing a time-reversal-invariant version of a parent model for fractional Chern insulators 43,44 .
Electronic interactions and spontaneous quantum anomalous Hall effect. The tunable realization of isolated time-reversal symmetric topological flat bands in twisted ZrS 2 is an ideal starting point for the stabilization of a host of correlated topological states of matter, ranging from interaction-induced quantum anomalous Hall effects 45 to elusive fractional Chern and topological insulators 46-49 . To investigate the role of electronic interactions and propensity for correlated topological phases in twisted ZrS 2 , we now study the topmost moiré Kagome band at fractional fillings and augment the effective three-band Kagome tight-binding description [ Fig. 3]-derived from the continuum theory, and a continuous function of the twist angle-via a screened Coulomb repulsion. The interaction is constrained for simplicity to a local Hubbard (U∑ ini"ni# ) and nearest-neighbor density (ðU 0 =2Þ∑ <ij>σσ 0n iσnjσ 0 ) interaction, expected to be a good approximation for screening due to metallic gates 50,51 . Suppose first that the top-most topological moiré Kagome band is tuned to half filling via electrostatic gating. A nontrivial spin Chern number precludes a straightforward Wannier tight-binding representation of this individual band. Instead, as deeper fully-filled valence bands are energetically separated, the low-energy behavior can be captured starting Kagome tight-binding model via projecting Coulomb interactions U; U 0 onto the Bloch states of a single fractionally-filled flat topological band, in direct analogy to lowest Landau level projections for the fractional quantum Hall effect. The resulting interacting problem is governed by an effective Hamilto-nianĤ whereĉ y kσ ;ĉ kσ create/annihilate electrons in the flat band with Bloch momenta k, ϵ k denotes the residual band dispersion, L is the system size, and Here, momenta k; k 0 ; q are defined in the moiré Brillouin zone, a i denote the moiré lattice vectors, and α, σ denote the sublattice and spin degrees of freedom. Since a sufficiently short-ranged interaction U > U 0 mainly imparts a local energetic penalty for electron pairs of opposite spin occupying the same Kagome "X" sites, a flat-band ferromagnetic instability generically ensues 52 at half filling of the top-most quantum spin Hall band, in direct analogy to quantum Hall ferromagnetism 53 . The resulting spontaneous spin-polarized state is gapped and aligned in the z direction-it exhibits a quantum anomalous Hall effect by virtue of filling a quantum spin Hall band for one spin component only; this fully spin-polarized state is an exact gapped zero-energy ground state in the absence of dispersion 45 and entails a quantized Hall conductivity. To study its robustness to the finite residual dispersion of the moiré band, we evaluate the phase diagram via exact diagonalization of Eq. (5) on a 4 × 4 unit cell cluster. Figure 4a depicts the phase diagram as a function of twist angle (which parameterizes the lowenergy electronic band structure and Bloch states) and interaction strength U vs bandwidth W of the top-most band. A robust quantum anomalous Hall state emerges for interactions on the order of four times the moiré bandwidth and remains robust over a wide range of twist angles. Notably, the underlying mechanism is distinct from the observed quantum anomalous Hall effect in twisted bilayer graphene, relying instead on the intrinsic topologically non-trivial moiré band structure due to strong spin-orbit coupling and obviating the necessity for concurrent valley polarization and substrate effects.
Fractional Chern insulators at fractional fillings. Having established a readily-accessible interaction-induced quantum anomalous Hall state at half filling, we now turn to the possibility of more exotic topologically-ordered phases at fractional fillings.
Here, previous theoretical works established numerous examples of topological tight-binding models 46,54-57 that a fractionally filled flat band with non-zero Chern number can in principle behave analogous to a fractionally filled Landau level and realize Abelian and non-Abelian fractional quantum Hall states in the absence of an external magnetic field. Uniformity of the Berry curvature is a key figure of merit 58 and determines, jointly with the ideal droplet condition for the Fubini-Study metric 59-61 , the propensity for flat Chern bands to host fractional quantum Hall phases. Figure 3f quantifies Berry curvature fluctuations ðΔΩÞ 2 ¼ R BZ ½ΩðkÞ À ffiffi ffi 3 p C=8π 2 2 dk for the Kagome moiré bands, where C ¼ ±1 is the spin Chern number. One finds that the Berry curvature flattens monotonically as the twist angle is reduced, with fluctuations substantially suppressed for the third Kagome valence band at small angles.
To demonstrate twisted ZrS 2 as a candidate to observe FQH physics with external magnetic fields, we focus on the conceptually simplest Laughlin ν = 1/3 state at 1/6 hole doping, and study the interacting problem at small twist angles in exact diagonalization. Analogous to the half-filled case, electrons in the almost-flat band can avoid local Coulomb repulsion U via spontaneous spin polarization, yielding a robust ferromagnetic instability as a function of U [ Fig. 4b, right axis, dashed line] over all investigated twist angles. However, spontaneous spin polarization due to U now leaves a single Chern band at 1/3 hole doping, with the resulting electronic phase governed by longerranged Coulomb interactions U 0 . To study the propensity to realize a Laughlin state, we numerically investigate the resulting phase diagram as function of bandwidth W=U 0 [Fig. 4b, left axis]. For W = 0, corresponding to the Landau level limit of a perfectlyflat Chern band, our exact diagonalization calculations for 6 × 5 unit cells reveal a three-fold ground state degeneracy for periodic boundary conditions [ Fig. 4c] with a gap to well-separated manybody excitations which persists as a function of system size. These ground states lie in three total momentum sectors that match the generalized Pauli principle for FCIs 56 , flow into each other upon adiabatic insertion of a magnetic flux through handles of the torus (periodic boundary conditions) and remain energetically separated from excitations, confirming the ν = 1/3 FCI 44,56 . Combined, these results indicate the robust stabilization of a fractional Chern insulator. The conclusions remain largely independent of the twist angle, and the FCI persists upon inclusion of finite band dispersion W until the many-body excitation gap closes for W=U 0 $ 0:3 [ Fig. 4b, false color].

Discussion
Having established a robust correlated quantum anomalous Hall phase at half filling and evidence for a ν = 1/3 fractional Chern insulator at one-sixth hole doping, an interesting follow-up question concerns the role of proximal deeper moiré valence bands, beyond the single-band approximation. For interactions that exceed the single-particle gap to other bands but remain smaller than the overall bandwidth of the three Kagome bands, the robustness of fractional Chern insulator phases has been welldocumented 62 , in direct analogy to Landau level mixing in the conventional quantum Hall effect. A more substantial challenge however stems from details of possible longer-ranged electron interactions and exchange processes, which could serve to either enhance or suppress the stability of the fractionalized phases at different filling fractions. These processes sensitively depend on the screening environment and gating 50 , and microscopic calculations present a substantial methodological obstacle for twisted materials 63,64 . Conversely, analyzing the potential stability of more exotic yet more fragile non-Abelian quantum Hall states remains an interesting topic for future investigation. Furthermore, for sufficiently small twist angles, if the Coulomb repulsion exceeds the overall bandwidth of the three Kagome bands, sufficient screening could serve to form a local moment at overall half filling ν = 3/2. Such a Kagome Mott insulator would constitute a Moiré realization of a paradigmatic frustrated magnetic model, which has been under intense scrutiny for the potential to host an elusive quantum spin liquid phase.
Beyond the (fractional) quantum anomalous Hall effect, the realization of flat-band quantum spin Hall insulators further opens up the possibility to realize a myriad of unconventional ordered states of matter with non-trivial topology, including time-reversal invariant fractionalized phases, or topological superconductors. Consequently, twisted ZrS 2 bilayers constitute a promising and tunable materials platform for such investigations, granting access to a novel and exotic regime of ultra-strong spin-orbit coupling that is not readily realizable in conventional crystalline solid-state systems. More broadly, a natural question concerns the extension of similar ideas of pseudospin potential engineering and strong spinorbit coupling to other transition-metal dichalcogenide heterostructures such as TiS 2 and HfS 2 with a multi-component character of the valence band edge. At the same time, the emergence of a moiré Kagome lattice from the fortuitous but robust interplay of geometry and interlayer coupling at small twist angles opens up a new pathway towards a moiré realization of magnetic phases in a paradigmatic frustrated system.

Methods
First-principles calculations. Ab initio calculations are performed with the Vienna Ab initio Simulation Package (VASP) 65 based on density functional theory (DFT). Plane-wave basis sets are employed with an energy cutoff of 450 eV. The pseudopotentials are constructed with the projector augmented wave method 66 and the exchange-correlation functionals are treated within the generalized gradient approximation (GGA) 67 . Only the Γ point is considered in the calculations due to the large size of the moiré supercells. A vacuum region larger than 15 Angstrom along the z-axis is applied to eliminate artificial interactions between periodic slab images. All atoms are relaxed until the forces on each atom are less than 0.01 eV/ Angstrom. Van der Waals corrections are applied with the Tkatchenko-Scheffler method 68 during the relaxation. The figures for the atomic structures and the charge density distributions are generated with the VESTA code 69 .
In this work, DFT calculations are only used to provide a reliable single-particle description of the top moiré valence bands of twisted bilayer ZrS 2 . The role of strong electronic interactions within the flat band is subsequently investigated in detail starting from the continuum theory and the effective tight-binding models described in the main text (with parameters extracted from the DFT calculations), using largescale many-body exact diagonalization calculations (see below). Although the band gaps in the systems to the conduction band are underestimated by the DFT calculations at the GGA level, the conduction band lies at high energies and remains empty, hence does not affect the emergent many-body state at small hole doping of the flat moiré valence bands. Only the band dispersion and the shape of the moiré bands are relevant in this work and these are well captured by GGA. As shown in refs. 38,70 , many-body corrections to GGA for 2D transitional metal dichalcogenides mainly appear as a rigid shift of the bands such that band gap is enlarged.
Exact-diagonalization calculations. Fractional Chern insulating phases and the spontaneous quantum anomalous Hall effect are studied using exact diagonalization calculations of the many-body ground state of the projected interaction Hamiltonian [Eq. (5)]. Calculations are performed for L 1 × L 2 unit cell clusters with periodic boundary conditions and discrete momenta k = n 1 b 1 /L 1 + n 2 b 2 /L 2 , with n i = 0, …, L i − 1 and reciprocal lattice vectors b 1;2 ¼ ½2π; ± 2π= ffiffi ffi 3 p > . Electron spin is explicitly included. Simulations of the interaction-induced quantum anomalous Hall effect at half filling are performed for 4 × 4 clusters. Results for ν = 1/3 fractional Chern insulators at 1/6 hole doping are obtained for 24-site (6 × 4) and 30-site (6 × 5) clusters.

Data availability
The raw data sets used for the presented analysis within the current study are available from the corresponding authors on reasonable request.

Code availability
Custom codes used in this work can be provided by the corresponding author on reasonable request. Ab initio calculations were performed with VASP (version 5.4.4).